GeodeticShiftFromWGS84 Subroutine

private subroutine GeodeticShiftFromWGS84(input, WGS84lat, WGS84lon, latOut, lonOut)

shifts geodetic coordinates relative to WGS84 to geodetic coordinates relative to a given local datum.

Arguments

Type IntentOptional Attributes Name
type(Coordinate), intent(in) :: input
real(kind=float), intent(in) :: WGS84lat
real(kind=float), intent(in) :: WGS84lon
real(kind=float), intent(out) :: latOut
real(kind=float), intent(out) :: lonOut

Variables

Type Visibility Attributes Name Initial
real(kind=float), public, parameter :: MOLODENSKY_MAX = 89.75*degToRad

Polar limit

real(kind=float), public :: WGS84_a

Semi-major axis of WGS84 ellipsoid in meters

real(kind=float), public :: WGS84_f

Flattening of WGS84 ellisoid

real(kind=float), public :: WGS84hgt
real(kind=float), public :: a

Semi-major axis of ellipsoid in meters
real(kind=float), public :: da

Difference in semi-major axes
real(kind=float), public :: df

Difference in flattening
real(kind=float), public :: dx
real(kind=float), public :: dy
real(kind=float), public :: dz
real(kind=float), public :: f

Flattening of ellipsoid
real(kind=float), public :: hgtOut

Source Code

SUBROUTINE GeodeticShiftFromWGS84 &
!
(input, WGS84lat, WGS84lon, latOut, lonOut )

IMPLICIT NONE

! Arguments with intent(in):
TYPE (Coordinate), INTENT(IN) :: input
REAL (KIND = float), INTENT(IN) :: WGS84lat, WGS84lon

! Arguments with intent(out):
REAL (KIND = float), INTENT(OUT) :: latOut
REAL (KIND = float), INTENT(OUT) :: lonOut

! Local variables:
REAL (KIND = float) :: hgtOut
REAL (KIND = float) :: WGS84_a !! Semi-major axis of WGS84 ellipsoid in meters 
REAL (KIND = float) :: WGS84_f !! Flattening of WGS84 ellisoid    
REAL (KIND = float) :: a !!Semi-major axis of ellipsoid in meters
REAL (KIND = float) :: f !!Flattening of ellipsoid
REAL (KIND = float) :: da !!Difference in semi-major axes
REAL (KIND = float) :: df !!Difference in flattening
REAL (KIND = float) :: dx
REAL (KIND = float) :: dy
REAL (KIND = float) :: dz 
REAL (KIND = float) :: WGS84hgt


! Local parameters:
REAL (KIND = float), PARAMETER :: MOLODENSKY_MAX  = 89.75 * degToRad !! Polar limit          

!------------end of declaration------------------------------------------------

!initialize
WGS84hgt = 0.
WGS84_a = 6378137.0 
WGS84_f = 1. / 298.257223563
a = input % system % ellipsoid % a
f = input % system % ellipsoid % f

!If input datum is WGS84, just copy to output
IF ( input % system % datum % code == WGS84 ) THEN
  latOut = WGS84lat
  lonOut = WGS84lon
  RETURN
END IF

!If latitude is within limits, apply Molodensky 
IF ( WGS84lat <= MOLODENSKY_MAX .AND. &
     WGS84lat >= - MOLODENSKY_MAX  ) THEN
     
   da = a - WGS84_a
   df = f - WGS84_f
   dx = - input % system % datum % dx
   dy = - input % system % datum % dy
   dz = - input % system % datum % dz
   CALL MolodenskyShift(WGS84_a, da, WGS84_f, df, dx, dy, dz, WGS84lat, WGS84lon, &
                         WGS84hgt, latOut, lonOut, hgtOut)
   
                         
   !                      Molodensky_Shift(WGS84_a, da, WGS84_f, df, dx, dy, dz,
   !                            WGS84_Lat, WGS84_Lon, WGS84_Hgt, Lat_out, Lon_out, Hgt_out)
     
ELSE !apply 3 steps method

  WRITE (*,*) 'Latitude is outside Molodensky limits'
  WRITE (*,*) '3-step datum transformation not yet implemented'
  STOP

END IF


END SUBROUTINE GeodeticShiftFromWGS84